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Abstract 

A new approximate method to calculate exchange-correlation contributions in the framework of 
first-principles tight-binding molecular dynamics methods has been developed. In the proposed 
scheme on-site (off-site) exchange-correlation matrix elements are expressed as a one-center (two- 
center) term plus a correction due to the rest of the atoms. The one-center (two-center) term 
is evaluated directly, while the correction is calculated using a variation of the Sankey-Niklewski 
approach generalized for arbitrary atomic-like basis sets. The proposed scheme for exchange- 
correlation part permits the accurate and computationally efficient calculation of corresponding 
tight-binding matrices and atomic forces for complex systems. We calculate bulk properties of se- 
lected transition (W,Pd), noble (Au) or simple (Al) metals, a semiconductor (Si) and the transition 
metal oxide Ti02 with the new method to demonstrate its flexibility and good accuracy. 
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I. INTRODUCTION 



The application of first-principles simulation techniques is becoming a research tool of 
increasing importance in materials science, condensed matter physics and chemistry, and 
molecular physics and chemistry. Most of these techniques are based on Density Functional 
Theory (DFT) 2] which creates an important simplification of the many-body quantum 
mechanical problem. Typically, DFT calculations are performed within the Kohn-Sham 
approach using the Local Density Approximation (LDA)^ or a Generalized Gradient 
Approximation (GGA) 4]. These total energy quantum mechanical methods can be used 
to calculate forces on atoms, and thus perform first-principles molecular dynamics (MD) 
simulations. Such simulations have been very successful in the description of a variety of 
properties of different materials. But, in spite of the important simplifications introduced 
by DFT and related approximations (e.g. LDA,GGA), they still require a huge amount 
of computational resources. This problem has severely limited the range of applications of 
these simulation techniques to situations with small numbers of atoms (~ 100-200) in the 
unit-cell, and short simulation times. 

Due to the computational limitations, first-principles simulation techniques have been 
mainly directed to the study of the energetics and electronic structure of diverse materials, 
surfaces and molecules. Typically, a good guess for the atomic structure is obtained before 
the calculation, and the first-principles method is then used to refine the geometry, obtain the 
electronic structure, and to compare the total energy of a few competing structures. These 
methods, however, have been very rarely applied to elucidate complex atomic structures, 
that require the exploration of an extensive phase space of possibilities, when no a priori 
answer, or approximate good guess, is already available. More importantly, the application 
of first-principles methods to investigate complex kinetic processes in materials (e.g. the 
atomic motion of atoms on a surface, kinetic pathways, molecular reactions, etc) is still very 
limited, due to the computational resources required for these calculations. 

It is clear that the usefulness of first-principles simulation techniques can be greatly 
extended if appropriate approximations are made, with the purpose of increasing the com- 
putational efficiency, with as little loss of accuracy as possible . This idea 
has prompted the development of first-principles tight-binding molecular dynamics (TBMD) 

methods HQ 

y.llf|. whose main characteristics are: (1) a real- space technique (i.e. no need 
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for super-cells or grids), (2) optimized atomic-like orbitals 

efficient, two-dimensional, tabulation-interpolation schemes 1. 5] to obtain flic effective TB 
Hamiltonian matrix elements as well as their derivatives to obtain the forces. 

The main advantage of such techniques is computational efficiency which makes them 
ideal first-principles exploratory tools. The use of first-principles TBMD methods as a 
exploratory tool can be complemented with more accurate calculations, if necessary; once 
stimulating results or ideas are obtained, final results can be refined by performing more 
accurate and time-consuming calculations (plane-waves DFT - e.g.^J - or even many-body 
- e.g. [3 - calculations). 

In this paper we report on new developments for the treatment of exchange-correlation 
contributions in first-principles TBMD methods, and their implementation in the FIREBALL 

nnn 

code [11, [15|, [16|| . The approximations used so far to handle these contributions are analyzed in 
Section II. In Section III we present our new approach to calculate the exchange-correlation 
contributions, and in Section IV we present some results for several materials to illustrate 
the performance of the new approach. 



II. AB INITIO TIGHT BINDING: FIREBALL 

Fireball Q, Q Q is a first-principles TBMD simulation technique based on a self- 



12lUi| functional. The energy functional is written 



consistent version of the Harris- Foulkes 

as: 



EtotW)] = E £ ™ - E ee [p{r)] + E xc \p{r)) - J p{f)V xc [p{^]d 3 r + E v 



ion— ion j 



where p(r) is the input density, which will be allowed to vary, and will be determined 
selfconsistently. The first term is a sum over occupied eigenstates, e n , of the effective one- 
electron Hamiltonian, 

(- l -V 2 + V[p\^4, n = e n i) n] (2) 

the potential V is the sum of the ionic potential, Vi on (r), (typically represented by a pseu- 
dopotential), a Hartree potential, and an exchange-correlation potential V xc 

V\p\ = v im {?) + J P ~P^r + VUp(r)]. (3) 
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In equation (JIJ E ee is an average electron-electron energy, 

E M = l -U^mr', (4) 



Eion-ion, is the ion- ion interaction energy 



Ei on -i on ^ -, (5) 

ij | fi'j | 

(Zj is the nuclear or pseudopotential charge of atom i at position Ri), and -Ee C [p] is the 
exchange-correlation energy First-principles MD simulations can be performed once the 
forces 

F t = (6) 

on each atom i are evaluated. 

The efficiency of calculations based on the Harris functional is associated with the possi- 
bility to choose p(r) in the above equations as a sum of atomic-like densities, Pi(r) : 

p(r) = ^2Pi{r)- ( 7 ) 

i 

In the Fireball method, confined atomic-like orbitals are used as a basis set for the de- 
termination of the occupied eigenvalues and eigenvectors of the one-electron Hamiltonian, 
Eq. (j2J). The fireball orbitals, introduced by Sankey and Niklewski (SN) are obtained 
by solving the atomic problem with the boundary condition that the atomic orbitals vanish 
outside and at a predetermined radius r c where ip(f)\ r=r . c = (see Fig. 1). An important 
advantage of the fireball basis set is that the Hamiltonian (Eqn. (j2J)) and the overlap matrix 
elements are quite sparse for large systems. The electron density p(r) is written in terms of 
the fireball orbitals (pn m (r) = M (r) {i is the atomic site, / represents the atomic subshell - 
e.g. 3s, 3p, 3d, etc., and m is the magnetic quantum number) 

p(r)=E^I^H 2 - (8) 

In this way four-center integrals are not required for the calculation of the Hartree terms, 
and all the two- and three center interactions are tabulated beforehand and placed in in- 
terpolation data-tables which are no larger than two-dimensional Hamiltonian matrix 
elements are evaluated by looking up the necessary information from the data-tables. 
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In practice, the atomic densities pi 

Piif) = ^,<lilm\<t>ilm{r)? ■ (9) 

in i 

are approximated to be spherically symmetric around each atomic site % (i.e. qu m = qu m ')- 
Self-consistency is achieved by imposing that the output orbital charges q° ut (obtained from 
;he occupied eigenvectors ip n of equation (J2J)) and input orbital charges coincide (see ref. 
3 and Q :or further details). 

The remaining difficulty is the efficient calculation of exchange-correlation interactions 
within a first-principles TB scheme. One possibility is to use non-standard DFT and intro- 
duce the exchange-correlation energy and potential as a function of the orbital occupancies 



2l|. In this paper, however, we opt for the more traditional approach in which 
exchange-correlation contributions are calculated as a functional of the electron density 
p(r). Within this line, two different methods have been previously proposed for the prac- 
tical calculation of exchange-correlation terms, using data-tables similar to those for the 
Hartree contributions. These two methods are: 



A. Sankey-Niklewski approximation 

The basic idea introduced by SN [1] is to write down the non-linear in p(r) exchange- 
correlation matrix elements in terms of matrix elements of p(r) . These later matrix elements 
are easily tabulated in data-tables no larger than two-dimensional, similar to those required 
for the Hartree terms. 

Consider the matrix elements < ^ | V^ c [p] 1^)^ > of the exchange-correlation potential. For 
each matrix element < M |\4 c [p]|0j, >, expand V xc [p] in a Taylor's series 

v xc [ P ] ~ v xc [p, w ] + v;' c [/vKp - /v) + - ( 10 ) 

around an appropriate "average density" p M „: 

_ < (PM4>» > fu> . 

< q>iAq>v > 

With this choice of p^ the second term in the expansion for < ^^jV^cfp] \(f> v > is zero, and 
the next term is minimized This yields 

< <P^\V xc [p]\<p u >=< p\V xc [p]\u >~ V xc [p^} < fjt\v > (12) 



SN [1] realized that corrections to this approximation are required for the case < >=< 
\x\v >= 0, and they devised a scheme to add those corrections specifically for a minimal sp 3 
basis set. 

The first step towards improved an improved exchange-correlation contributions in our 
TBMD scheme (see Sec. Ill) will be, precisely, to extend their ideas for a general atomic-like 
basis set. Also, it will be shown below that the SN average density approximation is not 
accurate for the exchange-correlation-energy on-site terms < p\e xc \p > for atoms with a 
significant valence electron density such as transition metals (e.g. Au, Ag, Pd, etc.); in Sec. 
Ill we will present a new scheme that corrects this problem. 



B. Horsfield approximation 

An alternative approach to deal with exchange-correlation terms within a first-principles 
TBMD method was proposed by Horsfield j^j], who introduced a many-center expansion 
based on Eqn. (jZJ. In this approach we can distinguish two cases jij] (i M is the atomic site 
corresponding to orbital fi and % v corresponds to orbital v\ 

(a) = i v = i 

< p\V xc [p\\v >^< n\V xc [Pi]\v > < fi\{V xc [pi + pj] - V xc [pi\)\v >, (13) 

(b) (v = i) ^ (i„ = j) 

< n\V xc [p\\v >=< p\V xc [pi + pj}\u > + < tA^Vxc\pi + Pj + Pk] - V xc [pi + pj]\v > .(14) 

Although practical experience has shown that this is an accurate approach in many cases, 
the on-site terms ( case (a) ) are not always well-approximated by the above equation, and 
additional numerical integrals are necessary for those term P, |(|. Another shortcoming 
of this approach is the fact that most of the computational time required to create the 
data-tables within this approximation is spent in the calculation of the exchange-correlation 
terms, reducing the computational efficiency 
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III. A NEW EXCHANGE-CORRELATION SCHEME FOR AB INITIO TIGHT- 
BINDING 



We now present a new approximation to calculate exchange-correlation contributions in 
a first-principles TBMD method. This new approximation is composed of two main parts. 
First, we generalize the SN approximation for an arbitrary atomic-like basis set. Second, we 
propose an improved approximation which combines the best features of both the SN and 
Horsfield approximations. 

To generalize the SN approach beyond sp 3 basis sets, we define average densities p^ using 
auxiliary, spherically symmetric orbitals (fu (defined below). The use of these spherically 
symmetric orbitals solves problems associated with the zero overlap < p\v >= cases; 
moreover, this approach for calculating exchange-correlation terms is consistent with the 
spherical approximation used in calculating Hartree terms. 

We define spherically symmetric orbitals tpu = ip^ for each atomic subshell corre- 
sponding to atomic-like orbitals 4>u m as follows. First, we consider the atomic-like orbitals 

(fiiim = uj u (r)Y lm (tt) (15) 

where wu(r) is the radial part of <pu m and Yi m (Q) the spherical harmonic associated with the 
angular part. Next, we define the spherically symmetric orbitals by 

ipu = uu(r)Y 00 (Q) (16) 

where Co is the positive root of (oj) 2 (see Fig. 1); a) is defined this way in order to avoid 
spurious cancellations in the importance sampling calculation of p^ (Eq. ITTj) which occur 
with certain intra-atomic cases (e.g. two different s-orbitals on the same atom). 

With these auxiliary orbitals we now define average densities for each matrix element 
(/!, u) as 

= < ^> b > (17) 

< <pp\<p v > 

This new definition for the average densities p^ v (using the auxiliary orbitals <p instead of 
the atomic orbitals <ft) solves all problems related to zero overlap (< (p^u >= 0), since now 
< (Pultfv 0; moreover, the use of auxiliary orbitals represents an improvement in the 
"importance sampling" calculation of p^ u for the non-zero overlap cases. Regions of positive 
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FIG. 1: Two radial functions corresponding to Fireball <f>\ and phi2, and auxiliary ip2 orbital. The 
auxiliary orbital (p% (not shown) coincides with <f)\. 

overlap are no longer "artificially" canceled by regions of negative overlap: both positive 
and negative overlap regions add- up now in this new definition of p^. 

Note that since the orbitals (fa are spherically symmetric, the same value of is obtained 
for all the matrix elements (fi, v) associated with a given pair of atomic subshells (i, I; i', I'). 
Also with our new definition of p^ u , we have in general < p\p\v >^ p^ < p\v >, and we 
keep the next term in the Taylor's expansion for 14 c [p]: 

< p\V xc [p}\u >~ V xc [p^} < p\v > +VJJP/0,] (< p\p\v > -p„ v < p\v (18) 

We refer to this approximation as the generalized SN approach (GSN approach). 

The average density approximation is not very accurate for certain matrix elements, 
particularly < p\e xc \p >. This is not a serious problem (as we demonstrate in Sec. IV) 
for the determination of the structural and/or electronic properties, since basically only the 
absolute value of the total energy is affected (i.e. it represents a rigid shift in the total 
energy). Nevertheless, the next step in our new approximation is to correct this inaccuracy. 
For this purpose, we use the best features of the SN and the Horsfield approximations. As in 
the Horsfield scheme, we distinguish two cases: (a) on-site (i^ = and (b) off-site matrix 
elements. 

(a) = i u = %. As a first step in our approximation, we simply add and subtract 
a contribution associated with the atomic density p« at site i, and write, formally, the 
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matrix element as a one-center contribution plus a correction, in similarity with the Horsfield 
approach: 

< p\V xc [p\\v >=< n\V xc [pi]\u > +(< p\V xc [p]W > ~ < p\V xc [pi]\is >V (19) 

The one-center (first) term is easily calculated and tabulated, and we use the GSN approach 
discussed above to evaluate the correction: 

< p\V xc [p}\u > ~ < p\V xc [pi]\u > 

+ V xc [p^ u ] < p\v > +V! cc [p^ u } ^< p\p\v > -p^ < p\v 

- V xc [pi] <p\v> -V£ c \pi] (< p\p t \v > <p\v>) (20) 

with 

Pi = 1 ( 21 ) 

(indices p, v have been omitted in pi, for clarity). 

(b) (in = i) ytz [i v = j). Proceeding in a similar manner as for the on-site matrix elements, 
we obtain for the off-site matrix elements: 



< p\V xc [p]\v > = < p\V xc [pi +Pj]\l> > +( < l*\Vxe\p]W > ~ < P\V xc [pi + Pj]\v > ) (22) 

- < p\V xc [pi + pj]\u > 

+ V xc [p^} < p\v > +VjJp ltt ,](< p\p\v > -p^ < p\v > 

- V xc [pij] < p\v > -V^piM < n\{pi + pj)\u > -p^ <n\v>\ (23) 



with 

= lrfft±ftfe (24 ) 
< > 

(indices p,v omitted for clarity). In Eqs. (jSPI and (J23J) p^ u , which includes all density 
contributions, is defined using Eq. (fTTj) . 



IV. RESULTS 

In this section we present results illustrating the performance of our new exchange- 
correlation scheme discussed in section III. Transition metals contain a significant valence 
electron density (the <i-electrons) , mixed with a free-electron-like density (the sp-bands), 
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and thus represent good test cases for the different exchange-correlation schemes. In this 
section we present results for some transition metals (Au, Pd, W), and the transition metal 
oxide Ti0 2 ; we also show some results for a typical sp 3 metal (Al) and a semiconductor (Si). 




lattice parameter [A] 



FIG. 2: Total energy for bulk Au as a function of the lattice parameter a for different XC- 
approximations presented in this work. The energy scale on the right corresponds to the GSN 
approximation. 'This work' line shows the result obtained from Eqs. 1201 and l23l The dashed line 
is the result obtained using Eqs. EH and 1141 Basis set: sp 3 d 5 fireball orbitals with cut-off radii 
R c (s) = 4.6 a.u., R c (p) = 5.2 a.u. and R c (d) = 4.1 a.u. 

As an example we plot energy vs. lattice constant for Au. Fig. |2] shows the total energy 
of bulk Au as a function of the lattice parameter as calculated with the Fireball code 
using the new exchange-correlation scheme proposed in section III (i.e. Eqs. (j!2(Jj) and 
(J22J)). We used a sp 3 d 5 basis set with fireball orbitals defined by the following cut-off radii: 
R c (s) = 4.6 a.u., R c (p) = 5.2 a.u. and R c (d) = 4.1 a.u. Also we show in the results for 
GSN (equation (fTH)) h and Horsfield (obtained using Eqs. (fl3|) and (jl4j) . i.e. without the 
additional numerical integral for the on-site terms). These results demonstrate how critical 
it is for the transition metals to have a good description of the on-site XC-contributions. 
Note the nearly rigid shift of the GSN result (scale on the right of Fig. |2J) by ~ 10 eV. 
As mentioned in Sees. II and III, this is related to the inaccuracy of the average-density 
approximation (Eq. fTS)) for calculating the on-site energy terms < fi\e xc \fi >: Au atoms 
contain a large electron density ( ~ 10 <i-electrons plus 1 s-electron) in the valence band. 
For comparison, this shift is only ~ 0.7 eV in the case of bulk Al. 



10 



TABLE I: Equilibrium lattice constants a and bulk moduli B for selected elements obtained using 
Eqs. EH] and 1221 (This work), and Eq. ^1 (GSN) for the exchange-correlation LDA contributions. 
We calculated these with sp 3 (Al,Si) or sp 3 d 5 (transition metals) basis sets of fireball orbitals with 
cut-off radii R c (in a.u.) as indicated. Also shown are PW-LDA and experimental values. 



Name 


Structure 


Basis set 


a (A) 


B(GPa) 


This work 


GSN 


PW-LDA 


Expt. 


This work 


GSN 


PW-LDA 


Expt. 


Au 


fee 


s4.6-p5.2-d4.1 


4.14 


4.21 


4.06 


4.07 


210 


197 


205 


171 


Pd 


fee 


s4.6-p5.0-d4.0 


3.96 


3.98 


3.85 


3.89 


215 


294 


220 


181 


W 


bec 


s4.7-p5.2-d4.5 


3.18 


3.17 


3.14 


3.17 


347 


320 


333 


310 


Si 


zbd 


s4.8-p5.4 


5.46 


5.50 


5.38 


5.43 


109 


98 


96 


100 


Al 


fee 


s5.3-p5.7 


4.04 


4.09 


3.97 


4.05 


93 


85 


84 


76 



Table I shows the calculated lattice parameter a and Bulk modulus B (obtained using a 
Murnaghan equation of state EOS) for Au, as well as for other transition metals (Pd, W), 
Al (a typical free-electron-like metal) and for Si (a typical semiconductor). These results 
have been obtained using either minimal sp 3 basis sets (Al,Si) or sp 3 d 5 basis sets (transition 
metals). The experimental values 22J, and the plane- waves LDA values (PW-LDA), are also 
presented in Table I. This Table shows that with the new approach to introduce exchange- 
correlation contributions the experimental lattice constants a are reproduced within ~ 2% 
while the bulk moduli are slightly overestimated by ~ 15%. The agreement is improved 
when comparing with the PW-LDA. Since the accuracy of first-principles TBMD methods 
is mainly related with the quality of the atomic-like basis set, improvements on the results 
presented in Table I are to be expected with a better choice for the basis set, either by 
improving the sp 3 or _sv 3 d 5 orbitals and/or adding new orbitals to the basis set (e.g. double 
basis sets, etc. 



Figures IHJ4J5I show the band-structures for the transition metals Au, W and Pd. The 
comparison with more accurate calculations (e.g. see j^J) shows that the band-structures 
are well-approximated within the present first-principles TB-scheme. In this case, there is 
practically no difference between the band-structures obtained with the GSN or the new 
approach (equations (j2 



Tetragonal rutile structure TiC>2 belongs to the space group P4i/mnm, containing 6 
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FIG. 3: Band structure of Au fee. LDA exchange-correlation terms are calculated using the new 
approximation discussed in Sec. III. Basis set: sp 3 d 5 fireball orbitals with cut-off radii R c (s) = 4.6 
a.u., R c {p) = 5.2 a.u. and R c (d) = 4.1 a.u..The dashed line represents the Fermi level. 
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FIG. 4: Band structure of W bee (see also Fig. |2J). Basis set: sp 3 d 5 fireball orbitals with cut-off 
radii R c (s) = 4.7 a.u., R c (p) = 5.2 a.u. and R c {d) = 4.5 a.u.. 



atoms per unit cell. The structural parameters for rutile structure Ti02 have been deter- 
mined to a high degree of accuracy from the neutron diffraction experiments performed by 
Burdett et al. 2?| . We have calculated the structural parameters and electronic bandstruc- 
ture for Ti02 in the rutile structure using the Fireball code and the exchange-correlation 
approach discussed in section III. For these calculations we have used a sp 3 basis for oxygen 
with cut-off radii R c {s) = 3.6 a.u. and R c {p) = 4.1 a.u., while for Ti a basis set of sp 3 d 5 
orbitals was used, with cut-off radii R c (s) = 6.3 a.u., R c (p)=6.0 a.u. and R c (d) = 5.7 a.u. 
The optimal structure is obtained by minimizing the total energy of the rutile {P&x/mnni) 
structures with respect to the lattice parameters a, c, and internal parameter u. We perform 
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FIG. 5: Band structure of Pd fee (see also Fig. EJ). Basis set: sp 3 d 5 fireball orbitals with cut-off 
radii R c (s) = 4.6 a.u., R c (p) = 5.0 a.u. and R c (d) = 4.0 a.u.. 



this minimization by a two-step procedure as outlined in Ref. [2£j . Table |H] summarizes the 
comparison of our results to the experimentally determined structural and elastic param- 
eters in Ti02; results from other theoretical works are also listed. This table shows that 
our results for the structural properties of TiO"2 in the rutile structure are within 1% of the 
experimental results of Burdett, et al. f27j. From the integrated EOS, we obtain a value for 
the bulk modulus B of 206 GPa which agrees well with the experimental value of 211 GPa 
[29I ] . In addition, our results agree well with the calculated results of others j3, Q • 

TABLE II: Theoretical results for structural and elastic parameters for Ti02 in the rutile struc- 
ture. Comparisons are made between our results and experimental results for the volume V, 
lattice parameters a, c, internal parameter u, and bulk modulus B; zero subscript represents the 
experimental results 



, c, m 

y 





V/V 


a/a 


c/c 


u/uq 


B (GPa) 


B [29J (GPa) 


Present work 


0.994 


0.997 


0.999 


0.994 


206 


211 


Other calculation [28] 


1.039 


1.013 


1.002 


1.001 


240 




Other calculation [30] 


1.021 


0.999 


1.002 


0.998 


209 





Using our theoretically predicted equilibrium lattice parameters, we have calculated the 
self-consistent electronic band structure for rutile Ti02 depicted in Fig. |U1 along the high- 
symmetry directions of the irreducible Brillouin zone. Table IIHI gives a summary of our 
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results in comparison to experiment and other calculations for the detailed features of the 
band structure. The upper valence band is composed of C>2p orbitals and has a width of 5.75 
eV. These results are in agreement with the experimental values of 5.50 eV |3l| . The lower 
band is 1.89 eV wide. Our results are consistent with other calculations 28J, Q. The 
calculated direct band gap at L of 3.05 eV is in agreement with the reported experimental gap 
of 3.06 eV (3^]. The LDA generally underestimates the experimental band gap for insulators 
and semiconductors, and the band gap obtained from ah initio plane-wave calculations for 
Ti0 2 is ~ 2.0 eV j28j. This underestimating effect of LDA is compensated in our results 
because we use a local orbital basis set. We also find an indirect band gap from L to M 
which is smaller than the direct band gap by 0.13 eV. 



TABLE III: Comparison of our present work to the experimentally determined electronic properties 
for Ti02 in the rutile structure. Definition of listed quantities are as follows: (1) Eg - D is the 
direct bandgap (L to L), (2) Eg - ID is the indirect bandgap (L to M), (3) Eyg is the upper 
valence bandwidth, and (4) Eq2 s is the Oxygen 2s state bandwidth. 



Structure 


Eg - D (eV) 


E (eV) g - ID 


E VB (eV) 


E 2s (eV) 


Present 


3.05 


2.92 


5.75 


1.89 


Experiments 


3.06 [22] 




5.50 [21] 




Others 


1.78 f30], 2.00 [22] 


3.00 [30], 2.00 [28J 


6.22 [2Q], 5.7 [22] 


1.94 [30], 1.80 [28J 



V. SUMMARY 

In summary, we have presented a new approach to calculate exchange-correlation contri- 
butions in first-principles TBMD methods. After a brief presentation of the basic theoretical 
foundations and practical motivation for these techniques, we have discussed the different 
approximations ( SN [j| and Horsfield j^j ) used so far to calculate exchange-correlation 
terms in these methods, using standard DFT (e.g. LDA). Then, in section III, we propose a 
new approach that corrects the main deficiencies of previous approximations in a practical 
manner (keeping always in mind computational efficiency). In this approach, on-site (off- 
site) exchange-correlation matrix elements are formally written as a one-center (two-center) 
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FIG. 6: Band structures for Ti02 in the rutile structure. The valence-band maximum is taken as 
the zero of energy. 

term plus a correction due to the rest of the atoms. The one-center (two-center) term is 
evaluated (and tabulated) directly, while the correction is calculated using the SN approach. 
For this purpose, a general (i.e. for arbitrary atomic-like basis set) version of the GSN 
approach has been also developed. 

The new scheme has been tested for some materials using the Fireball code and minimal 
sp 3 (for Al, Si and O) or sp 3 d 5 (Au, Pd, W, Ti) basis sets. The results, presented in section 
IV, show the good accuracy of the present first-principles TBMD approach as compared 
with experiment and other accurate calculations. 

Acknowledgments 

P.J. gratefully acknowledges financial support by the Spanish Ministerio de Educacion, 
Cultura y Deportes. This work has been supported by the DGI-MCyT (Spain) under con- 
tracts MAT200 1-0665 and MAT200 1-0 1534. JPL and HW greatly acknowledges financial 
support from DOE Grant No. DE-FG02-03ER46059 and from the Center for the Simulation 
of Accidental Fires and Explosions (C-SAFE at the University of Utah), funded by the De- 
partment of Energy, Lawrence Livermore National Laboratory, under subcontract B341493. 



15 



OFS thanks the NSF (DMR 99- 86706) for support. 



[1] O. F. Sankey and D. J. Niklewski, Phys. Rev. B 40, 3979 (1989). 
[2] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 
[3] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 
[4] A.D.Becke, Phys. Rev. A 38, 3098 (1988). 
[5] A.P. Horsheld, Phys. Rev. B 56, 6594 (1997). 

[6] A.P. Horsheld and A.M. Bratkovsky, J. Phys.: Condens. Matter 12 Rl (2000). 
[7] D. Porezag et al., Phys. Rev. B 51, 12947 (1995); T. Freuenheim et al. Phys. Rev. B 52, 11492 
(1995). 

[8] P. Ordejon, E. Artacho, J.M. Soler, Phys. Rev. B 53, 10441 (1996); D. Sanchez-Portal, P. 

Ordejon, E. Artacho and J.M. Soler, Int. J. Quant. Chem. 65, 453 (1997). 
[9] P. Pou et al., Int. J. Quant. Chem. 91, 151 (2003); R. Oszwaldowski et al., J. Phys.: Cond. 
Matter 15, S2665 (2003). 
[10] J. Ortega, Computational Materials Science 12, 192 (1998). 
[11] H. Eschrig and I. Bergert, Phys. Stat. Sol. B 90, 621 (1978). 

[12] J. Junquera, O. Paz, D. Sanchez-Portal and E. Artacho, Phys. Rev. B 64, 235111 (2001). 
[13] J. Ortega, R. Perez and F. Flores and A. Levy Yeyati, J. Phys.: Condens. Matter 12, L21 
(2000). 

[14] J. Ortega, F. Flores and A. Levy Yeyati, Phys. Rev. B 58, 4584 (1998). 

[15] A. A. Demkov, J. Ortega, O.F. Sankey and M.P. Grumbach, Phys. Rev. B 52, 1618 (1995). 

[16] J.P. Lewis et al, Phys. Rev. B 64, 195103 (2001). 

[17] J. Harris, Phys. Rev. B 31, 1770 (1985). 

[18] W. Foulkes and R. Haydock, Phys. Rev. B 39, 12520 (1989). 

[19] J.P. Lewis, J. Pikus, Th. E. Cheatham, E.B. Starikov, H. Wang, J. Tomfohr, and O.F. Sankey, 

Phys. Stat. Sol. 233, 90 (2002). 
[20] F.J. Garcia-Vidal et al., Phys. Rev. B 50, 10537 (1994). 
[21] P.Pou et al. Phys. Rev. B 62, 4309 (2000). 

[22] C. Kittel, Introduction to Solid State Physics, 6th ed. (Wiley, New York, 1986). 
[23] S.D. Kenny, A.P. Horsheld and Hideaki Fujitani, Phys. Rev. B 62, 4899 (2000). 

16 



[24] E. Anglada, J.M. Soler, J. Junquera and E. Artacho, Phys. Rev. B 66, 205101 (2002). 
[25] T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004). 

[26] D.A. Papaconstantopoulos, Handbook of the Band Structure of Elemental Solids, (Plenum 
Press, New York, 1986). 

[27] J. K. Burdett, T. Hughbanks, G. J. Miller, J. W. Richardson, and J. V. Smith, I. Am. Chem. 

Soc. 109, 3639(1987). 
[28] K. M. Glassford and J. R Chelikowsky, Phys. Rev. B 46, 1284 (1992). 

[29] T. Arlt, M. Bermejo, M. A. Blanco, L. Gerward, J. Z. Jiang, J. S. Olsen, and J. M. Recio, 

Phys. Rev. B 61, 14414 (2000). 
[30] S. D. Mo amd W. Y. Ching, Phys. Rev. B 51, 13023 (1995). 

[31] S. P. Kowalczyk, F. R. McFeely, L Ley, V. T. Gritsyna, and D. A. Shirley, Solid State Commun. 
23, 161 (1977). 

[32] J. Pascual and H. Mathieu, Phys. Rev. B 18, 5606 (1978). 

[33] K. Vos, J. Phys. C: Solid State Phys. 10, 3917 (1977). 

[34] R. V. Kasowski and R. H. Tait, Phys. Rev. B 20, 5168 (1979). 

99 



17 



